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ABSTRACT 


A capability for solving elasto-plastic plate bending problems is 
developed using assumptions consistent with Kirchhoff plate theory. 

Both bending and extensional modes of deformation are admitted with the 
two modes becoming coupled as yielding proceeds. The material work- 
hardens and, consistent with the fundamental theory of elasto-plasticity, 
loading is incremental and local unloading is permitted to occur. 
Equilibrium solutions are obtained numerically by determination of the 
stationary point of a functional which is analogous to the potential 
strain energy. The stationary value of the functional for each load 
increment is efficiently obtained through use of the conjugate gradient 
method (CGM) at modest computer storage requirements. The CGM was 
modified to take advantage of initial search direction vectors to provide 
possible reductions in computing time from one load increment to the next. 

This technique is applied to the problem of a large centrally through 
cracked plate subject to remote circular bending. Comparison is drawn 
between two cases of the bending problem. The first neglects the possi- 
bility of crack face interference with bending, and the second includes 
a kinematic prohibition against the crack face from passing through the 
symmetry plane. Results are reported which isolate the effects of elasto- 
plastic flow and crack closure. 



CHAPTER I 


INTRODUCTION 

1.1 General Comments 

The stress and strain fields in a plate containing a perfectly sharp 
crack subject to a variety of loading conditions has been the principal 
focus of a great deal of past and present research. The solution to the 
elastic extensional problem was given by Inglis [1]* in 1913. It was 
not until 1951 that Williams [2], using Kirchhoff fourth order plate 
theory, first published a solution to the problem of elastic bending of 
a centrally cracked plate. In 1961, Knowles and Wang [3] solved this 
same centrally cracked plate bending problem using sixth order Reissner 
plate theory. The principal deficiency with either the Williams or the 
Knowles and Wang solutions to the bending problem is their inability to 
account for the crack closure phenomena, i.e., both solutions model the 
crack face as a free surface either in the Kirchhoff or Reissner sense. 
This allows the material to overlap onto itself on the compressive side 
as bending and crack face rotation proceeds. 

The bending/extension interaction problem has been mainly attempted 
through superposition of elastic bending and extensional solutions. The 
work of Folias [4] and Wynn [5] are noted in this area. Specifically, 
the crack closure problem does not exist if the extensional load is large 
enough to open the crack sufficiently to preclude contact of the crack 


*Numbers in parentheses indicate references listed in the 
Bibliography. 
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face on the compressive side due to the bending load. Unfortunately, 
this crack face interaction does exist in many practical instances and 
its effect needs to be quantified. 

The extension/bending problem becomes even more insideous when 
material yielding is present. Under these conditions, not only is the 
bending/extension behavior coupled through interaction of the closed 
crack face, it is also coupled through material property considerations. 

Elasto-plastic solutions are available for the extensional case in 
the form of numerical solutions [{6), (7), and (8)] and closed form 
solutions in the case of special plasticity assumptions [(9), (10), and 

(11) ]. Unfortunately, the bending problem has not been attacked with 
such success. The principal reasons for this is the additional compli- 
cations of crack closure, free surface approximations for Kirchhoff 
boundary conditions, and the elastic core that persists in plate bending 
problems for the case of work hardening materials. Gonzalez and Brinson 

[12] have applied the Dugdale Strip model to the bending problem with 
some success but their solutions are limited both due to lack of crack 
closure considerations and not accounting for an elastic core. 

1 .2 Scope of Research 

Ideally, the analysis of a plate with a centrally placed through 
crack subject to circular bending of such a magnitude to cause large 
amounts of yielding would include both the effects of crack closure and 
coupling between bending and extension as yielding proceeds. A first 
step in that direction is provided in this work. 



The bending theory used in this study is based on classic Kirchhoff 
plate assumptions. Higher order theories for the bending analysis were 
deemed prohibitive when considered in conjunction with elasto-plastic 
material properties. The general elasto-plastic extension/bending 
theory is developed through use of a functional whose minimum corresponds 
to the equilibrium state for the elasto-plastic plate. Determination 
of this minimum through a numerical procedure constitutes the finite 
element approach used. Using this capability, the circular bending of 
an elasto-plastic plate with a centrally located crack is treated. The 
problem is essentially broken down into two parts: The first problem 

using the classic boundary conditions for a stress free surface (in the 
Kirchhoff sense) of the crack face (i.e., no closure) and the second 
which modifies these boundary conditions to reflect the physically 
realistic phenomenon 0 f crack closure. This is accomplished by adding to 
the functional, constraint equations that require the normal displacement 
of the crack face on the compressive surface to be zero. This couples in- 
plane and bending behavior through the constraint of the crack face and 
not by extending exterior plate boundaries to relieve the crack face 
interference. 

It should be pointed out that no corrections of the results are to 
be made concerning crack face warpage due to either three dimensional effect 
at the crack tip or to the wedging effect as the crack closes. Experi- 
mental evidence [13] indicates that these effects are important for crack 
length to plate thickness ratios of less than four, for which these 



results do not apply. Also, no geometric non-linearities are intro- 
duced either through the bending or stretching theories. However, since 
a circular bending field is applied to the plate, a developable surface 
is expected and the geometric non-1 inearities due to bending, should be 
small. 

1 . 3 Summary of Principal Findings 

Principal findings of this study clearly show that the mechanics 
of a centrally cracked plate subject to circular bending is strongly 
influenced both by the effects of crack closure and el asto-plastic flow. 
Results are given for in-plane and transverse deflections, transverse 
slopes, neutral axis shifts, stress and strain distributions, growth 
of yielded zones, and stress intensity factors for plates with and 
without closure effects accounted for. Comparison of these results 
between the two problems allow reduction of systemic errors and even 
though Kirchhoff plate theory was used, a good indication of the effects of 
elasto-plasticity and crack closure was obtained. 



CHAPTER II 


REVIEW OF LITERATURE 

Many theoretical and experimental studies have been performed 
concerning the problem of a stationary through crack in a plate subject 
to bending and/or extensional loading. Although the bulk of these 
studies deals with elastic behavior only, elastic-plastic results are 
available for the extensional problem mostly in the form of numerical 
determination of the stress state or asymptotic studies of crack tip 
singularities in plastic materials. There has been relatively little 
work reported on the elasto-plastic bending problem. 

2.1 Elastic Bending and Extension 

The early work by Williams [2,14] dealt with the bending problem by 
considering a wedge-shaped section and applying boundary conditions along 
radial edges of the section. Williams solved the biharmonic equation 
v^w = 0 of Kirchhoff plate theory by assuming a solution of the form 
w(r,o) = r A+ ^F(o,x). He then determined the functional form F(o,a) which 
solved the differential equations and obtained eigen solutions appropriate 
for various boundary conditions. Williams [15] then made the analogy 
between ~t'he plane extension v 4 f = 0 and bending v^w =. 0_where_'t l is the. 

Airy stress function and using the developed techniques, solved for various 
boundary condition combinations. In this paper, he suggested that the 
wedge-shaped section for the extensional problem be extended to 360° for 
the free-free boundary condition case. These conditions represent a 
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if 


semi-infinite V-shaped notch in an infinite plate. In [16], Williams 
published the results of the V-shaped notch or crack for the extensional 
case using the technique developed in [15]. 

In 1957, Irwin [17] presented the results of using the Westergard 
stress function [18] as applied to the through cracked plate. The 
singular stresses at the crack tip for the extensional case are 


xx 


V ( 2-rrr ) 


~r cos j (1 - sin j sin ^r) 


yy 
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(2-1) 


The significant points of this work are the 1A/F singularity , the 
ratio of the normal stresses at the crack tip (unity), and the angular 
variation of the stresses. 
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Williams [19] continued his eigen solution technique and in 1961 
published his results employing fourth-order classical bending theory 
with Kirchhoff boundary conditions. His singular stresses were modified 
by Sih, Paris, and Erdogan [20] to identify a bending stress and shear 
stress intensity factor. Using this later notation K , these singular 

D 

stresses are (see Figure 2-1 for plate geometry) 
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Sih, Paris, and Erdogan [20] went on to calculate Kg for various 
crack configurations using the method of complex variables as applied 
to planar problems by Muskhel ishvil i [21]. 

Williams [19] results showed that the ratio of the normal stresses 
at the crack tip in the plane of the crack are not unity as in Irwin's 
extensional problem and, although the (1/r) 2 singularity was the same, 
the angular variation was different between the extensional and bending 


case. 
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Knowles and Wang [3] used sixth order Reissner plate theory and 
showed that the functional form of the stress distribution is the same 
as in the extensional case. Moreover, the stress ratio is unity as in 
the extensional problem. The difference between Williams [19] and Knowles 
and Wang's results is due to the use of approximate Kirchhoff boundary 
conditions for the free surface, while the Reissner theory allows for 
the satisfaction of three separate boundary conditions, i.e., vanishing 
normal moment, twisting moment, and shear stress on the crack face. 

In 1969, Hartranft and Si h [22] applied Reissner plate theory to an 
infinite plate of finite thickness and showed that the bending stress 
intensity was a function of plate thickness. For vanishingly thin plates, 
these results reduced to those of Knowles and Wang [3]. 

2.2 Experimental Investigations 

In all of the theoretical studies mentioned above, the crack face 
was taken to be stress free in either the Kirchhoff or the Reissner sense. 
In reality, however, this is not the case; on the compression side of the 
plate, the crack physically closes on itself. There is ample experimental 
evidence that such is the case. For instance, Erdogan et al . [23] found 
that even in cases of some extension normal to the crack surface, the 
crack closes upon itself on the compression side of the plate. Smith 
and Smith [13,24], using a photoelastic technique, also found a strong 
influence of crack closure on the stress results and in fact found that 
crack closure could produce a non-conservative effect of as much as 40%. 
Smith and Smith [24] indicated a shift in the neutral axis that was 
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associated with the crack closure phenomenon. They compared their 
experimental results with the theoretical work of Hartranft and Sih [22] 
and found that as three-dimensional effects became prominent (i.e., 
crack length becomes small compared to plate thickness) crack closure 
effects dropped off. 

2.3 Elasto-Plastic Bending and Extensional Solutions 

A considerable amount of effort has been directed toward the 
extension of an elasto-plastic plate containing a stationary crack. 

There is much current work being done attempting to identify and make 
use of parameters which characterize near tip singularities of elasto- 
plastic solutions from both analytical and numerical approaches. 

Due to the complexity that non-linear material behavior adds to 
this problem, most of the analytic efforts reported in the literature 
have made use of various idealizations such as perfect plasticity or 
power law hardening materials. For example, for rigid-perfectly plastic 
materials. Rice [25] has shown that the shear strain exhibits a 1/r crack 
tip singularity. Also, both Hutchinson [9], and Rice and Rosengren [10] 
have found near tip solutions to power law hardening materials which 
have stress and strain type singularity of n/(l+n) and l/(l+n) where n 
is the hardening exponent. Much has been" learned from elasto-plastic - - 

finite element solutions presented by Lee and Kobayashi [8], Swedlow 
et al. [6], Miyamoto et al . [7], and others concerning the growth and 
shape of plastic yield zones, transition between elastic and plastic 
behavior, and the stress and strain behavior in the vicinity of the 
crack tip. 
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Experimental evidence [27] has shown that the detailed shape of the 
stress-strai n curve is of great importance when considering crack 
geometries. This has led many investigators to develop numerical capa- 
bilities to handle materials with general work hardening properties 
(i.e., [6], [7], and [8]). Swedlow, Yang, and Williams [6] using finite 
element techniques, indicate that for elasto-plasticity, the singularity 
of stresses at the crack tip decreases with load, while that for strain 
increases. Efforts (Rice and Tracey [26] and Levy, et al . [28]) have 
been made to develop elasto-plastic singularity elements with limited 
success. This is mostly due to the changing nature of the singularity 
as the plasticity develops, the dependence of the singularity on the details 
of stress-strain curve, and the difference between stress and strain 
singularities at the crack. 

Very little theoretical or experimental work has been done concerning 
the elasto-plastic bending of plates containing stationary cracks. 

Gonzalez and Brinson [12] have extended the Dugdale extensional model 
to the bending case. They performed an experimental study and found that 
the plastic zone varied with thickness in both tension and compression 
at the crack tip with the compressive zone continuing along the crack 
front on the compressi.ve side of the plate. They also found experimentally 
a shift in the neutral axis due to the effects of crack closure which 
were not included in their theoretical considerations. 
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2 . 4 Combined Bending and Extension 

Theoretical solutions to the elastic combined extension/bending 
problem for the most part have been obtained by superposition of local 
stress .fields due to the extensional load onto local stress fields due 
to the bending moment [5], [44]. These solutions do not account for the 
crack closure phenomena unless the extensional load is large enough to 
preclude contact of the crack face on the compressive surface due to 
bending. Adjusting the extensional load and bending moment such that 
the crack face is always open was suggested by Wynn [5] as a method of 
eliminating the crack closure problem. However, as Wynn's work pointed 
out, there are many practical instances of combined extension/bending 
where crack closure is of major importance. No elasto-plastic results 
are reported in the literature for the case of combined bending and 
extensional loading. 


V 



CHAPTER III 


DISCRETE ELEMENT FORMULATION 


As pointed out in Chapter II, theoretical solutions are not yet 
available for the problem of a through stationary crack in a plate 
subjected to bending and extension where the effects of crack closure 
are important. Experimental evidence indicates that analytical tech- 
niques are needed that can adequately reflect the shift in neutral axis 
due to crack closure and couple the resulting in-plane loadings to the 
transverse plate behavior. In the elasto-plastic problem, this is 
further complicated by the coupling between stretching and bending that 
occurs as yielding proceeds. 

3.1 Solution Technique 

Traditionally, for elastic behavior, minimum potential and comple- 
mentary energy theorems have been used either to derive equilibrium 
equations with complete boundary conditions or to determine approximate 
solutions to complicated boundary value problems. These theorems are 
also used in the development of general numerical capabilities for the 
solution of_ elastic structural problems, namely, the finite element 
method. Also, by discretizing the structure in question and treating 
the potential energy expression as a functional to be minimized, it is 
possible to cast the equilibrium problem as a mathematical programming 
problem. The minimization process may then be carried out by any 
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number of numerical algorithms developed from the mathematical 
programming point of view. For non-linear material behavior, minimum 
potential and complementary energy theorems are not valid. However, for 
incremental plasticity, rate formulations analogous to the minimum 
potential and complementary energy theorems can be obtained. Concep- 
tually these new functionals may be treated similarly to the potential 
and/or complementary energy with the interpretation being applied on an 
incremental rather than a total state basis. 

The problem addressed here is that of a flat plate that is loaded 
in such a way as to cause both extentional and bending modes of defor- 
mation at least locally. These displacements are considered to be small 
enough so that kinematic assumptions consistent with Kirchhoff plate 
theory can be made. Although the deflections are considered to be small, 
the resulting stress field is allowed to be of such a magnitude to cause 
considerable yielding of the material. The theory for this elasto-plastic 
material behavior relating increments in stress and strain to increments 
of applied load has been obtained by coupling the Prandtl-Reuss relations 
with Drucker's work-hardening hypothesis, [6] and [29]. 

Due to the complexity of the field equations and boundary condition 
-describing the stretching-bending elasto-plastic. plate behavior, only the 
simplest cases can be solved in closed form. However, because the problem 
is basically quasi-linear (see Swedlow [29]), it is possible to cast 
the boundary value problem as a set of linear sub-problems to be solved 



in succession. Each sub-problem is associated with a load increment 
and accumulated values of stress and strain. The solution to the 
complete elasto-plastic boundary value problem then consists of 
solving a succession of linear, anisotropic, inhomogeneous elasticity 
problems. 

The approach taken in this study is to discretize the structure 
and solve the resulting equilibrium problem of an assemblage of elements 
approximating the total plate structure. A functional is developed 
whose minimum corresponds to the equilibrium point of an elasto-plastic 
plate subjected to in-plane and transverse loadings. The assemblage of 
the discretized approximations to this functional is then numerically 
minimized for a succession of loads incrementally applied allowing the 
plasticity effects to build up with each increment. 

The plate element used in this work was constructed using two 
dimensional cubic Hermite interpolation functions. Bogner, Fox, and 
Schmit [30] used these displacement functions to develop an elastic 
plate bending element. The attractive feature of using cubic Hermite 
interpolation polynomials as displacement functions is that they have 
the property of continuity of displacement and slope at intraelement 
boundaries. Tong and Pi an [31] have shown that such an element satisfies 
kinematic admissibility for the entire plate while Birkhoff, Schultz, 
and Varga [32] have shown that the error bounds are proportional to the 
mesh size squared. The extension of this elastic element to the inelastic 
case was done using deformation theory for the case of pure bending and 
no stretching by Fox and Stanton [33]. 
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The discretized element representation used here incorporates the 
kinematic admissibility properties of the bicubic Hermite interpolation 
polynomials in addition to the incremental plasticity constitutive 
relations. The spatial variation of the material properties caused by 
the yielding is approximated by bilinear Hermite interpolation between 
additional nodes interior to the element that partition the element into 
four sections. The thickness integration is accomplished by a Lobotav 
quadrature of nine points. 

Using the discretized representation for the plate, a quadratic 
functional is developed whose minimum corresponds to the equilibrium 
point for the particular displacement-equilibrium problem. A modified 
version of Beckman's [34] conjugate gradient quadratic programming 
algorithm is used to minimize the functional and obtain the equilibrium 
solution. The modifications used take direct advantage of initial guess 
of both solution and direction of search. These modifications become 
extremely important as far as elasto-plastic solutions are concerned. 
Here, the solution from the previous load increment can be used for the 
initial guess and direction vector for the next increment. 

The primary advantage of the energy search technique to solve the 
equil ibrium problem is that of computer ^storage. By calculating the 
total potential energy as the summation of the energy contributions from 
each element, very large problems can be solved without requiring the 
assembly of the master stiffness matrix. Also, bandwidth becomes 
unimportant as zeros are taken advantage of regardless of their position 
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in the matrix. In plasticity solutions where direct advantage of initial 
guess solutions may be made, significant reduction in computation time 
may be realized. 

In this Chapter, the plate element stiffness matrix is derived 
and the conjugate gradient method developed. 

3.2 Constitutive Relations 

Since the mathematical theory of plasticity has received extensive 
treatment in the literature (e.g., [29]), there is little need to 
present a detailed derivation here. Rather, an outline of the steps and 
assumptions necessary to obtain the elasto-plastic flow rule used in 
subsequent calculations is presented. 

A flow rule for elasto-plastic deformation can be derived based on 
the following assumptions [29]: 

1. elastic and plastic strains are separable; 

2. quasi -thermodynamic postulate known as Drucker's hypothesis 
is imposed; and, 

3. the existence of a surface is presumed in stress space, 
known as a loading function, which is of the form 

<P ( J£ > J 3 ) “ K ^ 0 


(3-1) 
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where -j, is an arbitrary' function having the dimensions of stress whose 


arguments are 


J 

2 



S . . 
ij Ji 


J = - S S S . 
3 3 ij jk 


(3-2) 


Here , S.. = a.. - 1 a, . 6 . . is the deviator of the stress tensor o . . . 
-ij ij 3 kk ij ij 

The scalar k is a function of the plastic strain energy density 
given by 


W<P) = [a..dej^ 
J ij IJ 


(3-3) 


From Drucker's hypothesis, the plastic strain increment vector 
can be written as 

d e (P) = A - 11 _ (3-4) 

IJ 

here A is a constant of proportionality yet to be determined. Under 
these circumstances, the plastic strain increment is 

(d) y2 

- - de ij- - = . 2^7 ♦ij. +k* do k£ 

The scalar y 2 is given by 


2 - JL 

’ (p) tij+i 

eq 


(3-5) 
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3$ 

^ij ~ 30^ 


U l 

- eg 

= aAp) 


The total strain increment can be divided into recoverable and 
irrecoverable portions. This suggests the superposition of elastic and 
plastic strain increments. 


dc ij = dt i* ) + dE ij ) 


The elastic portion is described by Hooke's law. Combining (3-5) 
and (3-6), the flow rule takes the form 


2yds ij = da i j ' 1^ da kk 6 ij + Y da k* 


The inverse of (3-7) is 


v ^ e k?, 

■ d ‘« + i*T' - *7^ 
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Equation (3-8) can be divided by dt and by noting that the material 

derivative da../dt can be written as 3a. /3t = a., for a quasi- 
ij ij iJ 

static process, it is possible to write the flow rule independent of 
the time scale [6]. Hence the constitutive relations Equation (3-8) 
take the form 


a . . 

1J 


E 


i jk£ e k£ 


or 



F -1 • 
i jk£°k£ 


(3-9) 


(Note: is given in Appendix A for plane stress.) 

Equation (3-9) has been derived for a monotonically increasing 
equivalent stress, equivalent strain curve. 

3.3 Variational Formulation 

Force equilibrium for a static process is expressed by 


a . . . + F . = 

1J ,J 1 


(3-10) 


For a quasi-static process in which the displacements are small and 
F.. is both small enough not to cause yielding and independent of time, 
differentiation with respect to time yields 
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Defining dir as the first variation in the displacement field 
corresponding to a change in applied forces and integrating a ^ sik 
over the volume requires that 


"fjj 6,i i dv = 0 


(3-12) 


Equation (3-12) can be integrated by parts and the divergence theorem 
used to change the volume integration to one along the surface S. 


fa.. <5u • n . dS - a . . 6u • ^ dV = 0 

J ij i j J TJ 1 »J 

S„ v 


(3-13) 


Cauchy's stress equilibrium expression T n = c.-ru is used along the 

I I J J 

surface to produce the equation 


o . . du . . dv 

ij i ,J 


\ T r .' Au. dS = 0 

J i l 


(3-14) 


The strain rate, displacement rate equations £..=%( u. . + u, .) can 

l j i » J J , i 

now be used resulting in a rate form analogous to the Principle of 
Virtual Work. 



6e. . dv 


f T n 

J l 


<Su. dS 

l 


(3-15) 


V 
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Using the constitutive Equations (3-9), the stationary principle 
required to determine equilibrium can be established. That is. 


J Wk^ij dV ■■ I T ? 5 “i ds (3 - 16) 

V. s u 

Defining the functional 

’ * -ij dV - I T ? "1 dS < 3 - 17 > 


and since E.^ and 
Equation (3-16) gives 


remain constant at an instant in time, 
rise to the stationary principle 


= 0 (3-18) 

The statement for Equation (3-18) is as follows: 

Theorem I . 

Of all displacement rates u.. satisfying the given boundary conditions, 
those which satisfy the equations of equilibrium are distinguished 
by a stationary (extreme) value of the functional n. 

For the case of elasticity, it is clear from Equation (3-17) that 
the functional tt reduces to the potential energy and that Theorem I is 
the minimum potential energy theorem. In the same sense that Equation (3-18) 
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is the analogue to the minimum potential energy theorem, a functional 
tt* can be developed to be analogous to the complementary energy theorem. 
This functional is 

"* '■ hi* - J ds (3 - ,9 > 

v S u 

Since the volume is fixed and C... and u. remain constant at an 

ijla i 

instant in time, the adjunct to Theorem I is 

<$TT* = 0 

Therefore, the second theorem is: 

Theorem II. 

Of all the stress rate tensor fields o.. that satisfy the equations 

* J 

of equilibrium and boundary conditions where traction rates are 
prescribed, the "actual" one is distinguished by a stationary (extreme) 
value of the functional 7r*. 

3.4 Discrete Element Representation 

The basis for the displacement solution is the determination of 
the stationary point associated with the functional in Equation (3-18). 
Viewing this functional as the summation of contributions from an 
assemblage of discrete elements, it is clear that the displacement state 
that minimizes Equation (3-18) corresponds to the equilibrium solution. 
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By assuming a displacement rate state over each discrete element, the 
volume integration can be completed with only the coefficients of the 
assumed displacement state as unknowns. Therefore, the most important 
phase in the development of the discrete element is the selection of 
the individual element displacement rate functions. As pointed out by 
Bogner [30] and Shi eh [36], the displacement rate functions must meet 
several criteria to guarantee that the stationary point of the functional 
of the assemblage of elements corresponds to the global stationary 
point for the complete structure. The formulation for the functional 
indicates the requirements that must be met by a displacement rate 
function. These requirements are: 

1. Satisfy global kinematic admissibility by satisfying boundary 
conditions and enforcing intraelement continuity of displace- 
ment and slope rates. 

2. It must have at least as many undetermined parameters as 
degrees of freedom. 

3. The assumed displacement rate shape should be complete in 

the sense that increasing the number of elements monotonical ly 
increases accuracy of results. 

4. It must -allow- for rigid body displacements.otherwise element 
equilibrium could be violated. 

As indicated by Bogner [30], Stanton and Schmit [37], and others, the 
use of displacement rate states written as the sum of two-point Hermite 
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interpolation polynomials make the satisfaction of kinematic admissibility 
very convenient. By allowing for two in-plane and one out-of-plane 
displacement rates, there are 12 constants and 12 degrees of freedom per 
corner of a rectangular element. Using the two-point bicubic Hermite 
interpolation polynomials, these degrees of freedom are matched exactly 
with those of the polynomials. Also, rigid body displacements are 
accounted for exactly. As Key [38] has shown for the special case of 
pure bending of a flat plate, the element displacement rates are complete 
and assure monotonic convergence. 

For these reasons, the displacements rate shapes considered here 
are written as two-point Hermite interpolation polynomials. 

When values of f(x) and df(x)/dx = f'(x) are known at two points, 
then, using Hermite interpolation, 

2 

f(x) = ^ [H 0i (x)f(x i } + H li (x)f/(x i )] (3-21) 

i=l 

The one dimensional cubic Hermite interpolation polynomials 
H Qi (x), H^.(x) may be found in Hildebrand [39]. For an interval 
Os x s a, they are 

Hqi (x) = 1/a 3 [2x 3 - 3ax 2 + a 3 ] 

Hi i (x) = 1/a 2 [x 3 - 2ax 2 + a 2 x] 

H 02 (x) = - 1/a 3 [2x 3 - 3ax 2 ] 

H-^U) = Va 2 [x 3 - ax 2 ] 


(3-22) 
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For two-dimensional interpolates, products of two, one-dimensional 
interpol ation formulas are formed to describe the function over the 
rectangular domain of interest. The undetermined parameters are the 
values of the function, first derivatives with respect to the arguments, 
and a cross derivative. 


3.5 Rectangular Bending-Membrane Plate Element 

For the plate analysis to follow, several assumptions are made in 
addition to those already mentioned. These are: 

1. Geometrically, deflections are small compared to plate 
dimensions. 

2. The plate is thin and a state of plane stress exists, i.e., 
° zz * t X 2 = T yz * 0 (see Figure 3-1). 

3. Lines orthogonal to the middle surface of the plate before 
bending, remain so after bending. 

Expanding Equation (3-17) for the case of plane stress consistent 
with assumption 2 yields 


a b h/2 

* ' i J j % IE 11 * 2E 12 Vy + E 22 V + 2E 13 £ x 
o o -h/2 


• 2 


(3-23) 


2E 23 x xy 'y + E 33 O dxd * dz ' I Tj «Cl. dS 
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Taking advantage of assumption 3, the strain-rate displacement-rate 
equations are 



Using bicubic Hermite interpolation formulas for the displacement 
shapes, the displacement rate state for the rectangular element over 
region Osx s a, O^ysb where the unknowns are corner parameters of 
a typical element is (see Figure 3-1 for numbering scheme of i,j). 

I 

au. . 

u(x,y) = u. . H (x) H (y) + _li H (x) H (y) 
ij Oi Oj 8x^1 OJ 


I 


2 - 

9U . . 9 U . . 

+ __IiH (x) H (y) + — -liH (x) H (y) 
3y Oi Ij 9X9y 1i lj 


(3-25) 


This equation can be rewritten in matrix form as 

u(x,y) = (u} T {P(x,y)l 


(3-26) 


Complete vectors { u } and (P(x,y)} are given in Appendix B. 
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Combining (3-23), (3-24), and (3-26) results in 


r . T 
u 


* * 


w 


' k UU 

I 

i 

1 

, uv 
k 

k uw 

_k u _v 

r 

i 

i 

_i_ 

> 

> 

k VW 

k uw 

i 

i 

, vw 
k 

, ww 
k 

L 

i 




C ~\ r T 

u j 


- 




v >{B n } 


w 

V J 


(3-27) 


where { B n } is the vector of external forces applied on the nodes of 
the element, and each of the k 1J are 16 x 16 symmetric entries defined 
by the volume integrals as follows: 


CO ■ jl E n p x p x + 2E 13 p x p J + E 33 p y p y 


v - 


dv 


[k uv ] = 


| 2£ 12 P x P J + 2E 13 P x p I + 2E 23 t P y P P + 2E 33 P y P x 


dv 


[k uw ] = 


J I 2ZE 11 P..P T . + 2ZE-, oP..pJ... + 4ZE 1 - 5 P V P T 


V 


'1 2 x yy 


11 X XX 
1 3 r y r xx 

[ k "J -_A\hz r A* 2E 23 P x P y * E 33 P x P x 


+ 2ZE„P..pL + 2ZE 23 P y P yy + 2ZE 33 P y P iy 


13 x r X y 

dv 


dv 


[k vw ] - - | 


v 


2ZE 1 2 P y P xx + 2ze 1 3 p x p xx + 2ZE 22 p y p yy 


+ 2ZE 23 P x P yy + 4ZE 23 P y P xy + 2ZE 33 P x P xy 


dv 
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[k WW ] = 


Teh z 2 p p t 

I II xx x 


xx + 2E 12 z2p xx P Jy 


+ 4E 13 z2p xx P xy + E 22 z2p yy P yy 

+ 4E 23 Z Vxy + "SjVly 


-dv 


In these expressions the P's are the Hermite polynomial vectors 
listed in Appendix B. The subscripts on P indicate differentiation with 
respect to particular variables. 

Since the material properties are in general functions of 
both the stress state and spatial coordinates, special consideration of 
them is due before the volume integration can be performed. To 
accomplish the x-y integration, the element is broken up into four 
subregions. Values of 

h/2 

E i j(x,y) = j E i j (x ,y)dZ 

-h/2 
h/2 

C^(x,y) = J ZE..(x,y)dZ (3-29) 

-h/2 
h/2 

D-J^z.y) = J Z 2 E.^.(x,y)dZ 
-h/2 

k s region number as per Figure 3-1 
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are subject to bilinear Hermite interpolation within each subregion 
where the values of the four adjacent corners are the unknown parameters. 
This permits closed form integration in the x-y plane with numerical 
evaluation of E^j, C^, required only at the nine connecting 
nodes of the four subregions. This numerical integration in the thickness 
direction is accomplished using the Lobotav quadrature with nine thickness 
points. For the elastic case, k uw = k vW = [0]. 

Equation (3-27) represents the contribution to tt for the n t * 1 
element. The total value of the functional tt is the sum of u 11 over 
all the elements or 

NEL 

V -1 n 

n = ^ tt where NEL = number of elements (3-30) 

i = l 

The equilibrium position for the plate is the stationary value of the 
functional v. 

3.6 Minimization Algorithm 

The total value of the functional u is obtained by summing over all 
elements their respective contributions to it as indicated in Equation 
(3-27) . By cal l ing the element degress of freedom {-Y} and element .. . _ 
stiffness matrices [K n ], Equations (3-27) and (3-30) can be rewritten as 


= % { Y n } T [Kj{ Y n } - { Y n } T { D n } 


(3-31) 
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Defining a transformation of the element degrees of freedom {Y n } to 
the system degrees of freedom {X} such that 


{Y n } = [a n ]{X> 


(3-32) 


it is possible to write Equation (3-30) in the form 


NEL 

= £ TT n = v? {X} T [A]{X} - {X} T {F} 

n=l 


where 


and 


NEL 

[A] > £ [c* n ] T [K n ][a n J 

n=l 


{F} = 


NEL 

L 

n=l 


£ [a n ] T {B n } 


(3-33) 


Since [A] is a positive definite symmetric matrix, the stationary point 
of the quadratic- functional in (3-33) is a minimum. _ It _is now possible 
to cast the equilibrium problem in terms of a mathematical programming 
problem, i.e.. 

Minimize: *(x) = % {X} T [A]{X} - {X} T {F} 

Subject: No constraints. 


(3-34) 
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In this representation, {X} is the state variable vector of dimension 
N (number of system degrees of freedom), [A] is a positive definite 
N x N matrix, and { F } is the N-dimensional work equivalent force vector. 

There are many techniques for minimizing Equation (3-34). Essen- 
tially, they fall into two categories: direct search methods and 

gradient methods. Most success in structural mechanics problems has been 
through the use of gradient techniques [30], [35], [37], [40]. Although 
Gisvold and Moe [40] applied Powell's (see [42]) method of direct search 
in solving buckling problems with satisfactory results, it has been 
pointed out by both Box [42] and Beveridge and Schechter [41] that 
gradient methods rather than direct search methods are best suited for 
unconstrained minimum problems. This is especially true when an analytic 
form of the gradient of the objective function can be determined. It 
is for these reasons that a modified version of the conjugate gradient 
method as outlined in Beckman's paper [34] is used in this work. 

All gradient techniques locate the solution {HI to the minimization 
problem as the limit of a sequence { X Q } , { X-j } , { Xg where { X q } is 
an initial approximation to { H> . (Note that subscripts refer to particular 
vectors and not their components.) The movement from one approximation 
{ X /> to the next approximation {X-j+y}- proceeds along, some, specif ied_ 
direction {p^} a distance a.. That is 


{X i+1 l = (X 1 } + ai { Pi } 


(3-35) 
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The techniques of determining and {p^} comprise the primary 
differences between the various methods. By constructing the {p^} to 
be linearly independent, one forms a basis for the solution space { H } 
and hence, in linear problems, convergence to the solution in N steps 
is guaranteed. Beckman [34] constructs the { p^ } ' s by a Gram-Schmidt 
orthogonal ization process using the gradient of tt evaluated at the {X^l's. 
Hence the name, conjugate gradient. The distance of travel is 
determined by minimizing tt(X^ + -j) with respect to . This gives the 
maximum distance of travel along any direction vector. 

It should be noted that in this method as in the conjugate gradient 
methods used by Stanton and Schmit [37],.Bogner [3 3, Bogner, et al . [30], 
and Fox and Stanton [33], the initial direction vector is in the direction 
of steepest descents. A common characteristic of this starting direction 
is that the number of iterations to convergence for a particular problem 
is independent of the initial guess { X Q ) and hence { p Q > . This behavior 
was found by Beckman [34] and Stanton and Schmit [37]. Although this 
characteristic may be an advantage in many instances, for incremental 
elasto-plastic problems, a distinct advantage could be gained by making 
use of the solution at the previous load step as a starting point for the 
present load step. 

The following method is based on Beckman's method but instead of 
starting in steepest descent direction, an arbitrary initial direction 
is used. Constructing the solution vector sequentially as 

(H> = (X o > + a 0 lp 0 > + f Pi >' + “ 2 tP2 } + ••• + a N_i ^ Pn-1 ^ 
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it is possible to size a . such that the distance a., in the direction 
(p.j } minimizes the objective function n. This is if 

<W = { V + W 

Then is such that 


d '< X Ul> 

dot • 
i 


= 0 


Carrying out this operation gives 

(P i } T {R i > 

011 {p-j } T [A]{pi } 


where the residual vector { } is defined as 

(R i } = { F} - [A]{X.J 

It is now necessary to pick a new direction vector. This vector is 
determined by an [A] - orthogonal ization of { p Q } , {R-|}, •••> t ^n-1 ^ 

using a Gram-Schmidt process. This results in a set of linearly inde- 
pendent vectors spanning the N-space describing the solution vector of 
Equation (3-34). 
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The algorithm described above may be summarized as follows: 

1. {X Q } = An arbitrary initial guess to {HI. 

{R 0 1 = (BJ - [A]{X q > 

{ p Q } = an arbitrary initial guess to the direction from 

{X Q > to (X^, or {R }. 

= <Pj> T <V 
{p,) T tA](P 1 ) 

(X 1+1 ) - {*,) + 

3. (R i+] } = {F} - [A]{X i+1 } 

or 

{R i+] } = (R i ) - o i [A]{p i } 

{R i+1 } T [A]{Pi} 

ZJ. t £ = — - 1 ■ " - — — ■■ ■■ 

1 (P i ) T [A]{p i } 

5- (P i+1 ) - {R i+1 } + 3 i {P i > 

Steps 2 thru 5 are now repeated until the solution to the desired 
accuracy is obtained. This takes at most N-steps for linear problems. 
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The choosing of the initial guess {X Q > and direction { p Q > is 
very critical in this routine. However, for linear problems, choosing 
{ X Q } to be {0} and { p Q } to be proportional to the solution, results in 
the exact solution in only one iteration. This is easily seen as a 

o 

then simply scales { p Q > to the correct magnitude. For elasto-plastic 
analysis, this is very helpful for it has been found that from one load 
step to the next, the displacement states are roughly proportional. 

Hence, using the old solution for both { X Q } and { p Q } should pjovide very 
rapid convergence behavior. For elastic problems, approximate solutions 
are often available that may be used as { X Q > and { p Q } . If not, { p Q } may 
either remain arbitrary or determined by steepest descents with convergence 
still guaranteed in at most N-iterations . 

In the course of their work using the conjugate gradient method for 
geometric non-linear structural problems, Fox and Stanton [33] uncovered 
a strong scaling effect significant in the convergence of their problems. 
Essentially, the difficulty stemmed from large variation in diagonal terms 
of the [A] matrix. They found that by making all diagonal terms equal 
(say 1.) and the sum of off-diagonal terms as small as possible, convergence 
of the conjugate gradient algorithm was greatly increased. This behavior 
was suggested by Gerschgorin ' s disk theorem. Geometrical ly , this scaling 
transforms the N-dimensional hyperellipsoid described by tt ( X ) = constant 
into a more hyperspheroidal shape. The scaling transformation is 
accomplished by premultiplying the state vector {X} by a diagonal matrix 
whose diagonal entries are d. ^ = l./VA~j\ 
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Convergence may be determined by one or both of two methods. 

Either the magnitude of the residual vector may be forced to be very 
small or convergence to a certain digit accuracy of the objective function 
may be required. In practice, both methods can work well. The minimum 
residual vector criteria seems to be more restrictive, however, its 
value does not necessarily decrease monotonical ly from iteration to 
iteration. On the other hand, the objective function does converge 
monotonically to the minimum value and is thus attractive as a convergence 
criteria. 

It should be pointed out that linear constraint equations of the 

form 

[C]{X) + {D} = 0 (3-41) 

can be introduced into Equation (3-34) through the use of Lagrange 
multiplier concepts. For the linear constraint equations of Equation 
(3-41), Equation (3-34) can be written as 

tt = %{X} T [A]{X} - {X} T {B} - {A} T [C]{X} - {A} T {D}. (3-42) 

or 

= %{X} T [A]{X} - {X} T {B} 


TT 


(3-43) 
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where 



For m linear constraint equations* [C] is a nxm connection matrix, {A} 
is an m dimensional Lagrange multiplier vector, and { D} is an m 
dimensional constraint vector where n is the dimension of the original 
{ x } vector. With zeroes along the main diagonal, [A] is no longer positive 
definite. However, [A]^[A] is positive definite and a new conjugate 
gradient method for this case can be derived (see Beckman [34]) 
following the same reasoning as above. 

3.7 Numerical Implementation 

The procedures outlined in the previous sections were assembled into 
a computer program called PLATER [47]. An outline of the program organi- 
zation is shown in Figure 3-2. As input to the program, material properties, 
element geometry, and intraelement communication data are required. 

Boundary conditions and constraint equations are necessary on an incre- 
mental basis. The material properties include elastic constants and 
discrete point values of octahedral stress vs. octahedral plastic strain. 
Figure 3-3 shows the PLATER element along with nodal load and displacement 
sign conventions. 

From an operational standpoint, the program evaluates the constitutive 
relations of Equation (3-7) for each element at the beginning of each 
load step. The properties are then considered constant over the subsequent 



FIGURE 3-2 


PLATER Program Flow Diagram 










FIGURE 3-3 

PLATER NODAL LOADS AND DISPLACEMENT 
SIGN CONVENTION 
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note: 

POSITIVE NODAL LOADS ARE 
IN THE SAME DIRECTION AS THE 
CORRESPONDING DISPLACEMENT 
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loading increment. Elastic unloading is systematically accounted for in 
the procedure. Isotropic hardening is assumed for subsequent yielding. 

The program was verified by comparing available closed form solutions 
with numerical results from PLATER. A variety of boundary conditions 
were used for bending, extensional , and combined bending/extension 
problems. Table 3-1 lists the verification problems and gives some 
results for stresses and displacements. An indication of the discreti- 
zation error is given in Table 3-1 by comparing the results as the 
number of elements is increased. Obviously, certain types of problems 
require higher degrees of discretization than others. For example, 
resolution of stress results under point loads is poor until the number 
of elements in the vicinity of the load gets large. However, as 
Table 3-1 indicates, as the number of elements in a given area was 
increased, the results (displacements and stresses) converged monotonical ly 
toward the correct answer. 

The convergence properties of the conjugate gradient method (CGM) 
were also investigated extensively. Particular attention was paid both 
to ill-conditioned master stiffness matrices and to the effects of 
initial guess values of both {x Q }andip o }. In this context, ill-conditioned 
matrices imply those which have large weighting of off diagonal terms. 
Ill-conditioned matrices were examined since it was expected that zeroes 
would appear on the main diagonal and that there would be large weighting 
of off diagonal terms due to extension/bending stiffness representations. 
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TABLE 3-1 

Verification Solutions for PLATER 


ELASTIC PROBLEMS 


Bending Results 



Circular Bending 
Square Plate 

Simply Supported Square Plate 

Clamped Square 
Plate 




Center 

Load 

Uniform Load 

Center Load 

Number 

of 

W 1 

Jmax 

M x| 

W Jmax 

M x| 

w Jmax 

M x| 

w Jmax 

El ements 


J max 

M 

0 

Pa 2 , 0 3 

D 

J max 

10 2 

Pa 2 

2 , 
aa_ 10 3 

J max 

, o 2 

2 

qa 

t £'° 3 

1 

,0 

1 .0 

11.11 


HI 

1 

5.32 

2 

1.0 

1.0 

Ewai 

— 


4.98 

5.37 

4 

1.0 

1.0 

D El 

1.65 

uSSm 

4.80 

5.47 

8 

1.0 

1 .0 

BBS 

1 .88 

4.03 

4.79 

5.54 

EXACT* 

1.0 

1 .0 


1 .88 

4.06 

4.78 

5.60 

Solution 



m 






Number 

of 


Extension/Bending Results 


Circular Bending Plus Uniform Compression 


Elements 

^max 

max 

^max 


°yj 





-'membrane 

J BENDING 

1 

+.12 x 10" 2 

-.36 x 10" 3 

.0605 

1000. 

3000. 

EXACT 

+.12 x 10" 2 

-.36 x 10‘ 3 

.0605 

1000. 

3000. 


* NAVIER Solution, see Ref. [43]. 
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TABLE 3-1 (con't) 
Verification Solutions for PLATER 


ELASTO- PLASTIC PROBLEMS 


Bilinear Stress-Strain Curve; E^ = 10 x 10^, E^^ = 40 x 10^ 
Circular Bending || Applied Twisting Edge Moment 



PLATER 

EXACT 


PLATER 

EXACT 

^ield 

0 

5490 

5490 

^yield 

0 

3130 

3130 

wyield 

max 

1 .82 

1 .82 

^ield 

max 

1.95 

1 .95 


The Following Data is for 
M o * 8125 T 


The Following Data is for 
Edge Twisting Moment of M° y = 4630 


WJ 


max 


4.85 


4.85 


wj 


max 


5.55 


5.55 
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Results of these studies indicate the following: 

1. Convergence of the CG!i is strongly influenced by ill- 
conditioning of the [A] matrix. This can be somewhat 
alleviated by the scaling transformation suggested by Fox 
and Stanton [3 3]. 

2. The initial solution guess vector {X Q } has very little effect 
on the convergence properties. 

3. The selection of the initial search direction vector { p Q > 

has a strong influence on convergence. However, the selection 
of { p Q } for a multidimensional problem is very difficult and the 
steepest descents direction (used in standard CGil formulations) 
is better than a bad guess for { p Q } . Much advantage can be 
obtained where good approximate {p Q } vectors can be found from 
previous solutions. In particular, incremental plasticity 
solutions can take direct advantage of this property by 
constructing { p 0 } from one load increment to the next. 

4. CGH convergence properties can be greatly influenced by judicious 
choice in orientation and numbering of both system degrees of 
freedom and element descriptions. The more uniform the element 
map, the better the convergence. No attempt to optimize a 
particular element map was made. Non-uniformity of the element 
map has a much greater effect on the convergence properties than 
having elements with large aspect ratios involved. 



CHAPTER IV 


STATIONARY THROUGH CRACK IN A PLATE 
SUBJECT TO CIRCULAR BENDING 


The elasto-plastic plate bending capability, PLATER, is now applied 
to the study of the through cracked plate. The plate is to be of sufficient 
dimension as to model an infinite plate where the crack length is large 
compared to the plate thickness. The plate width and length to crack 
length ratios are to be large enough so that the plate exterior boundaries 
do not interfere with the crack behavior. Essentially, two problems 
are to be analyzed. The first one neglects the effects of crack closure 
and the second includes these effects. Comparison of the two solutions 
permits elimination of any systemic errors resulting from the numerical 
procedures and isolates the closure phenomena. Geometrically, both 
problems are dimensioned as show in Figure 4-1. The boundary conditions 
for the two problems are identified in this chapter. 

The exterior dimensions are set so that the finite plate models as 
closely as possible an infinite plate. No corrections for finite plate 
effects are made. Also, no corrections are made for crack face warping. 
Although crack face warping may be expected physically, the Kirchhoff 
assumption of planes remaining plane preclude the realization of the 
warping effect. However, as indicated by Smith and Smith [13], at crack 
length to plate thickness ratios greater than four, three dimensional effects 
become less important compared to crack closure effects. As shown in 
Figure 4-1, the plate under consideration does have 2a/h = 4.0. 
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FIGURE 4-1 

CRACKED PLATE GEOMETRY 


47 



0>P 



48 


The material properties are listed in Figure 4-2. The work hardening 

material chosen is felt to be rather typical of many common materials. 

The Ramsberg-Osgood power law is cited only for convenient comparison 

to other published results [ 6 ] as the PLATER program uses discrete 

points on the t . vs. curve. This curve is shown in Figure 4-2. 

oct 'oct a 


4.1 Boundary Conditions Without Closure 

The boundary conditions for the first problem neglecting the crack 
closure behavior are identical with the ones used by Williams [16], i.e., 
the crack face is modeled as a free surface while the exterior plate 
boundaries are subjected to constant applied normal moment. 

The free surface boundary conditions in Kirchhoff plate theory 
are obtained by coalescing three conditions (i.e., normal moment, twisting 
moment, and average shear force equal to zero) into two conditions 
admissible by fourth order differential theory [43]. These conditions 
along a line parallel to they axis are 

(M x ) • 0 

x=a 


V 

x 


3M XV 

<Q X - * 0 

x 9 y x=a 


(4-1) 


Also, reaction forces at the crack tip and plate corners warrent particular 
attention in Kirchhoff theory. As shown in Appendix D, no additional 
contributions to the boundary conditions are required for this analysis. 



FIGURE 4-2 
MATERIAL PROPERTIES 
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The complete boundary conditions for the plate with no closure 
considerations are as follows: 

Along Y axis. 


X = 0, Y = Y 

N = 0 

M 

= M 


X 

X 


N = 0 

V 

= 0 


xy 

X 


Along X axis, 

X = X, Y = 0 


N = 0 

y 




Along line X = A, 
X = A, Y = Y 
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3U , elV 

ay ax 


2z 
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axay 


0 


Along line Y = B 
0 <; X ^ A-a, Y = B 




au. , dy_ ? aw 
ay ax “ z axay 


o 


Along crack face. 
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The applied moment, l'i 0 , is applied incrementally by scaling the 
first elastic load so that the most sensitive element is brought to the 
proportional limit. The next three load steps are at 5% of this initial 
yield load. The remaining load steps are applied at 10% the initial 
yield load. Table 4-1 lists the load increment and accumulated load 
level for the no closure problem. 

4.2 Boundary Conditions With Closure Effects 

The problem described in the above section is to be reanalyzed 
including the effects of the crack face interference as the crack face 
rotates due to applied moments. Clearly, the only boundary conditions to 
change are those along the crack face. For this problem, the crack 
is no longer traction free but subject to a normal compressive load due to 

i 

the crack face interaction which prevents material overlap as the 
plate is bent. Since the PLATER program is quasi-three dimensional (the 
through thickness behavior being assumed in the Kirchhoff sense) the 
wedging effects of the compressive surface must be modeled in terms of 
mid-surface quantities. The constraint of the crack face on the compressive 
surface can be obtained from the Kirchhoff displacement relations of 
Equation (3-24). Referring to Figure 4-3, the displacement in the Y- 
direction (v) of the compressive contact edge (Z = - h/2) is required to 
be zero. The Kirchhoff displacement relation is 

= „ (x ,y ) - Z 

S ay 


v(x,y) 


(4-3) 



52 


TABLE 4-1 


Incremental Load Steps and Accumulated 
Values for the No Closure Problem 
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Along the crack face, 

v(x,B) = 0 = v $ (x,B) + | ^ (x,B) (4-4) 

Therefore the constraint equations are 

{ V + f{f}'° (4 - 5 > 

When these constraint equations are used in PLATER, the in-plane exten- 
sion and transverse bending behavior are coupled. The constraint equations 
are introduced into PLATER through the use of Lagrange multiplier concepts. 
The resulting Lagrange multipliers are given the interpretation of in- 
plane loads and transverse moments. Although the mid-surface line along 
the crack face is actually free, these Lagrange multiplier couples give 
rise to pseudo-stresses in the results. Physically, this result is 
very realistic with the neutral surface shifting reflecting the wedging 
effect of the crack face. 

Loading was incremental for this problem as it was for the no 
closure problem. Table 4-2 lists the load increments and accumulated 
values for the closure problem. 



TABLE 4-2 


Incremental Load Steps and Accumulated 
Load Values for the Closure Problem 



Applied Edge Moment 

Accumul ated 



n IN-LBS 

Val ue 

| gfl 1 ■ 


IN 


IN-LBS 


m 



L IN J 


i 

370. 

370. 

2 

18.5 

389. 

3 

18.5 

407. 

4 

37 

444. 

5 

37 

481. 

6 

37 

518. 

7 

37 

555. 

8 

37 

592. 

9 

37 

629. 

10 

37 

666. 

11 

37 

703. 

12 

37 

740. 

13 

37 

777. 

14 

37 

814 . 

15 

37 

851. 

16 

37 

888. 

17 

37 

925. 

18 

37 

962. 

19 

37 

999. 

20 

37 

1036. 

21 

37 

1073. 

22 

37 

1110. 

23 

37 

1147. 

24 

37 

1184. 

25 

37 

1221. 
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4 . 3 Discretization of the Crack Problems 

Geometrically, both problems have 77 elements and 96 nodes with 
the elements in the vicinity of the crack tip being (a/16) square. For 
the no closure problem where bending variables only are considered, there 
are 353 system degrees of freedom. For the closure problem, both bending 
and extensional modes of deformation are excited which results in 1046 
degrees of freedom plus 9 constraint equations. 



CHAPTER V 


NUMERICAL RESULTS 

The numerical data extracted from the PLATER analysis of the 
centrally cracked plate bending problems are put into graphical form for 
ease in distinguishing the effects of crack closure and material yielding. 
The detailed elastic results reported are from the first load increment 
of the elasto-plastic incremental solutions normalized to a common load 
level. The elastic results for the no closure case allow not only for 
direct comparison with the closure case but for detailed verification of 
the finite element solution with the analytic solutions of Williams [19], 
Equation (2-2). 

5.1 Elastic Displacement Data 

The deformed shapes of the elastic plate with and without closure 
effects are shown in Figure (5-1). These figures clearly indicate that 
the kinematic boundary conditions on the plate were met. Of particular 
interest is that the in-plane displacement normal to the crack face along 
the compressive surface of the plate was essentially zero (on the order of 
10 in.). Also of importance is that the crack face remains flat, 

consistent with Kirchhoff plate theory. 

Figures (5-2) are contour plots of the transverse displacement, w. 
These plots have been normalized to the quantity M Q /2D(l+v) which is the 
coefficient of the (x +y^) solution to the circular bending problem of a 
plate without a crack. These contours are circles whose radius is the 


FIGURE 5-1 

ELASTICALLY DEFORMED PLATE 
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FIGURE 5-2 

CONTOURS OF TRANSVERSE DEFLECTION W [2D ( I +u)/M 0 ], 

FOR ELASTIC PROBLEMS 
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square root of the quantity plotted. As can be seen, the presence of the 
crack changes the circle into an ellipse around the crack. The ellipse 
for the crack closure model has a larger minor axis than that for no 
closure. Note that the transverse displacement is not a particularly 
sensitive parameter in studying the closure phenomena and that the 
elliptical contours of the transverse displacement caused by the crack, 
quickly (within three plate thicknesses) degenerate into circular contours. 
This circular contour is the developable surface of a sphere which is 
generated by the circular bending of a plate without a crack. Hence, it 
is expected that the errors due to neglecting large transverse displace- 
ments are small, especially since the transverse displacements are small 
(compared to the plate thickness) in the vicinity of the crack tip. In 
all cases, the slopes are small (approximately 0.03). 

Figures (5-3) and (5-4) are contour plots of the transverse slopes ^ 
and . These quantities are normalized to M Q /D(l+v) so that comparison 
with the solution of the plate with no crack can be made. Note that 
without a crack, is straight line parallel to the Y axis and ^ is a 

bX dj 

straight line parallel to the X axis. Deviations from this straight line 
behavior are indicative of the crack presence in the analysis. Comparisons 
of Figures (5-4a) and (5-4b) shows that ^ is affected by the crack 
presence but not markedly different with and without the closure effects 
present. Since this is the slope parallel to the crack face, this result 
is expected. On the other hand, Figures (5- 4a) and (5-4b) show that the 
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FIGURE 5-4 

CONTOURS OF TRANSVERSE SLOPE dw/dy [D ( I +v )/M 0 ], 

FOR ELASTIC PROBLEMS 
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transverse slope normal to the crack face is strongly influenced by the 

presence of the crack and exhibits a very pronounced change due to crack 

closure effects. Hence, is a very sensitive measure of the crack 

ay 

closure phenomenon. Again, and indicate that the region of influence 

oX ay 

of the crack is about 3h. 

The in-plane displacement normal to the crack face is plotted in 
Figure (5-5). A plan view of the tensile side of the plate is shown so 
that the overlap of material is indicated by negative displacements. This 
is the case for the no closure model with the compressive edge being 
overlapped onto itself (shown in Figure (S-5)) and the centroidal surface 
being at zero displacement. The closure model shows the compressive 
edge being at zero displacement and the centroidal surface being displaced 
exactly one-half of the value on the tension surface. The crack opening 
displacement is seen to be higher for the closure model. This implies a 
higher stress intensity factor if a linear elastic fracture mechanics 
approach is taken. (For further discussion of this, see Section 5.3.) 

5.2 Elastic Stress Results 

In Figures (5-6a), (5-7a), and (5-8a), the moment components, (M , 

A 

My, M^y) are plotted vs. r/a along the e = 0 axis for the elastic response. 

These moment distributions for the no closure problem compare favorably 

to the Williams [19] solution near the crack tip. For example, M x and 
-k 

exhibit the (r) 2 singularity with M x /M going negative near the crack 

y 

tip as a consequence of Kirchhoff boundary conditions. Also present is 



CRACK OPENING DISPLACEMENT (IN.) 
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FIGURE 5-5 

ELASTIC CRACK OPENING 
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FIGURE 5-6 

ELASTIC MOMENT (M x /M 0 ) AND 
STRESS DISTRIBUTIONS (cr x /cr 0 ) 65 



r/a 

a) MOMENT DISTRIBUTION ALONG 6*0° AXIS 
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FIGURE 5-7 

ELASTIC MOMENT (M y /M 0 ) AND STRESS 
(<r y /o- 0 ) DISTRIBUTIONS 
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FIGURE 5-8 67 

TWISTING MOMENT (M xy /M 0 ) AND SHEAR STRESS 
(r xy /<r 0 ) DISTRIBUTIONS FOR ELASTIC RESPONSE 



a) TWISTING MOMENT DISTRIBUTION ALONG 0*0° AXIS 
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the singular behavior of M vu at the crack tip as approached along the 
crack face. This is a consequence of the approximate nature of the free 
surface assumption in Kirchhoff theory. It should be noted that the force 
boundary conditions remote from the crack are satisfied and that the zero 
normal moment along the crack face is also satisfied for the no closure 
problem. 

Figure (5-7a) shows the non-zero normal moment along the crack 
face for the closure problem. This moment is the result of the wedging 
of the crack face during crack closure. As shown in Figure (5-9), 
there is actually an in-plane load normal to the crack face generated 
as a result of the crack face interference along the compressive edge 
during bending. Transferring this load to the centroidal surface of the 
plate (Figure (5-9)) requires the addition of a couple to satisfy equili- 
brium. As shown in Figure (5-9), although the mid-surface is actually 
a free surface, there are these pseudo forces and moments calculated at 
the mid-plane of the plate. The shift in the neutral axis resulting 
from the in-plane load is further discussed in Section 5.4. 

In Figures (5-6b,c,d), (5-7b,c,d) and (5-8b,c,d), the stress distri- 
butions through the thickness are shown for (r/a) = 3/32 and e = 0°, 90°, 
180°. For the no closure case, these plots show the symmetrical stress 
distribution (about the mid-plane) required by pure bending. Also, the 
linearity of stress and the resolution of the zero a at the crack face 
for no closure (e = 180°) is noted. The shift in neutral surface and the 
moment at the crack face is evident in the closure case. 
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FIGURE 5-9 

VARIATION OF IN-PLANE MEMBRANE LOADS 
FOR THE ELASTIC CLOSURE PROBLEM 

9 
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The variation of the stress components with the angle around the 
crack tip are compared with the Williams solution in Figures (5-10), 

(5-11), and (5-12). Fairly good agreement is achieved with the usual 
difficulty in resolving finite element stress results across inter- 
element boundaries being present. 

5.3 Fracture Mechanics Interpretation 

As discussed in Section 5.1, the elastic results may be given a 
linear elastic fracture mechanics interpretation. This approach requires, 
among other things, the extraction of the bending stress intensity data 
from the finite element results. Since displacement data are more accurate 
than stress data in finite element studies, the displacement results are 
compared with the Williams 09] solution in terms of the bending stress 
intensity factor. For convenience, data along the crack face are used. 

That is, using the expression from 09] for the transverse deflection w 
results in in-plane displacements u = -z|^ and v = -z— of the form 


- & b^Ccos e(- cos + 3 lj~ v ) . cos ~) - sin e(sin sin |] 


1 * 3/0 
2zb 2 r[cos e(cos2e + ^-) + sin e sin 2 q] - z Y r 


(5-1) 


- | zb^r^sin e[- cos ~ cos |] + cos e[sin ~ [y^ y sin f]> 


1 v 3 J 2 

- 2zrb 2 sin e[cos 2e + y^-] - cos e sin 2e - zar ' 



STRESS RATIO 
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FIGURE 5-10 

STRESS VARIATION WITH ANGLE COMPARING THE NO 
CLOSURE SOLUTION WITH WILLIAMS [l 9] SOLUTION 




MOMENT RATIO M v /M 
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FIGURE 5-11 

MOMENT VARIATION WITH ANGLE 
FOR ELASTIC RESPONSE 



FIGURE 5-12 

MOMENT DISTRIBUTION ALONG 6*90° 
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Along the crack face where e = 180°, the expressions for u and v at 
z ~ - h/2 are 


2hh 2 r 3/ 2 

u = — + yr 

1+v 


= + 


6hb-j r 
7+v 


1/2 


+ a r 


3/2 


where b 2 , y,a are constants and 


b 


1 


2(74-v)(1+v) k 

3(3+v)EhVi7 B 


K 


B 



(5-2) 


Dividing both sides of Equations (5-2) by V r/a, it is possible to get 
the v displacement in the form 


where 



(5-3) 


a is a constant coefficient for the higher order terms of the 
Williams solution, and 

= 40+v) K 

E(3+v)V2tt B 


e 
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As in all finite element studies, the results near the crack tip are not 

as good as those away from the crack tip. However, by plotting v/Vr/a 

from the finite element results and extrapolating to the crack tip, 

values of a and 3 are obtained from which Kg is calculated. Figure (5-1 3a) 

shows results of this procedure for the no closure model with the value 

2 k 

obtained from the extrapolation being 14% lower than Kg = 6M Q /h (ira) -5 . 

This is indicative of the systemic error of the finite element model and 
is unimportant since a comparison between the closure/no closure model is 
of prime interest. Figure (5-1 3b) shows the closure results obtained 
through the same process. Figure (5-14) compares the results through the 
thickness for the closure/no closure problems. Note that K g is zero at 
the compressive edge for the closure problem and at the mid-plane for the 
no closure problem. The value of Kg is also negative on the compressive 
edge in the no closure problem. This leads to ambiguity as to the ptysical 
interpretation if any, that can be applied to the bending stress intensity. 

The linearity is due to the Kirchhoff assumptions. The Kg is about 20% 
higher on the tension surface for the closure problem than for the no 
closure case. This is trend is fairly consistent with experimental observation 
in photoelastic materials [24], Actually, from photoelastic studies. 

Smith and Smith [24] obtained 40% higher values of Kg but they could not 
entirely isolate the closure effects in their experiments because of a 
number of other events that operate simultaneously (e.g., three- 
dimensional effects such as crack extension and crack face warping.) 
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FIGURE 5-13 76 

DETERMINATION OF ELASTIC BENDING STRESS 
INTENSITY FACTORS ON THE TENSION SURFACE 
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FIGURE 5-14 

COMPARISON OF BENDING STRESS INTENSITY DATA 
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5.4 Neutral Surface Shifts 

The strongest indicator of the crack closure phenomenon is the shift 
in the neutral surface. Actually, the definition of a neutral surface 
is arbitrary in the sense that either zero forces, strains, or displace- 
ments may be used to determine it. Only kinematic quantities will be 
used here; specifically either e = 0 or v = 0 will define the neutral 
surfaces. 

In his work on the combined bending/extension elastic problem, 

Wynn [5] applied an extensional load at the exterior plate boundaries in 
such a magnitude as to preclude contact of the crack face surfaces 
during bending. Under these circumstances, crack closure is obviously 
no problem. However, adjusting the superimposed extensional load in the 
exact amount to relieve the contact on the compressive surface does not 
correctly model the crack closure phenomenon. This , process puts the plate 
in tension while the actual mechanics along the crack face applies a 
compressive loading along the line of contact. Also, these forces are 
variable along the crack face as the wedging actions of the crack 
closure changes. This difference in mechanism is indicated by the shift 
in neutral surface defined by = 0 shown in Figure (5-15). If an 
exterior tension had been applied, the neutral surface would shift toward 
the compressive side of the plate. However, the constraint against 
material overlap is guaranteed by contact of the compressive surface which 
adds a compressive force and a normal bending moment along the crack face. 
As shown in Figure (5-15), this results in a shift of neutral surface (as 
defined by = 0) toward the tension surface. 



FIGURE 5-15 

NEUTRAL SURFACE AS DEFINED BY 
=0 FOR THE ELASTIC CLOSURE PROBLEM 
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NEUTRAL SURFACE 
€ y = du/dy - 
Z d 2 w/dy 2 = 0 






80 


By defining the neutral surface as v = 0, the shift is toward the 
compressive surface as shown in Figure (5-16). In fact, the neutral 
surface is identical to the compressive contact line at the crack face, 
and decays toward the mid-plane away from the crack face. The neutral 
surface is coincident with the mid-plane at a distance of three plate 
thicknesses from the crack face. The maximum shift is at the plate center 
and gradually becomes less pronounced as the crack tip is approached. 
Figure (5-16) shows the neutral surface shift at various stations of x. 

5.5 Elasto-Plastic Results 

The growth of plastic yield zones is the most dramatic feature of 
the elasto-plastic results. Yielding is defined at those points for 
which tq C j vs. yqct 1S nonlinear. For the no closure problem, the yield 
zones are symmetric about the mid-plane as required by the Kirchhoff 
bending theory (see Figure (5-17)). Also a consequence of the Kirchhoff 
theory is the dominance of the singular shear stress term along the crack 
face approaching the crack tip. This dominates the plasticity behavior 
and forces the yield zones to generate along the crack face first and 
then to propogate into the material. Without the singular shear stress, 
it is presumed that the yield zone shape at any surface parallel to the 
mid-plane, would be similar to that found in extensions! plane. stress 
(see [6]). 

Comparing the no closure results with the analyses of Brinson and 
Gonzalez [12] and Brinson, et al . [48], using the Dugdale model, indicates 
that the elastic core present in this analysis is a very important feature 
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FIGURE 5-17 

YIELD ZONES FOR THE NO CLOSURE PROBLEMS 
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in the elasto-plastic bending of the plate. Being a rigid/perfectly 
plastic model, the Dugdale (or strip) model requires full penetration of 
the yield zones, i.e., no elastic core present at the crack tip. 

Figure (5-1 7b) shows that the yield zone shape through the thickness 
obtained from this analysis is not cusp shaped as in the strip model and 
that the plate has a definite elastic core. 

The formation of yield zones is completely different in the case 
of the closure model. As shown in Figure (5-18), the yield zones form on 
the tension surface first and proceed through the thickness before 
yielding on the compressive edge. It is noted in Figure (5-18a) that 
the in-plane yield zone growth appears to be a combination of the plane 
stress extensional case and the bending case with no closure. The 
important feature, though, is that yielding progresses on the tension 
surface first prior to initial yield on the compressive contact surface. 
The yield zones tend to grow along the crack face for the closure 
problem due to the compressive stress field built up by the contact of 
the crack face and the shear stress due to Kirchhoff boundary conditions. 
Although the stresses are a bit higher in the closure case, the yield 
zones are smaller. This is due to the greater constraint of the crack 
surface in the closure case causing compressive fields along the crack 
face which inhibits yielding. 

The through thickness stress distributions in Figure (5-19) show 
neutral axis shifts and nonlinear stress behavior as loading builds up. 
The crack closure results in shifting the neutral surface and in causing 
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FIGURE 5-19 

THROUGH THICKNESS STRESS DISTRIBUTIONS 
FOR THE ELASTO- PLASTIC PROBLEMS 
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initial yielding to occur on the tension surface. As yielding proceeds, 
the elastic core takes up more of the load until yielding is initiated 
on the compressive surface. 

The moment redistribution as loading proceeds is indicated in 
Figure (5-20). By comparing the closure case with the no closure case, it 
is seen that for a given load, more redistribution takes place with the 
closure model. This is consistent with the fact that for a given applied 
moment, less material is yielded due to the constraint of the crack face 
for the closure model. 

The neutral axis shift, as defined by the surface where v = 0, is not 
greatly affected by the yielding process. The neutral surface does shift 
back toward the mid-plane of the plate as the plasticity builds up, but 
this behavior is only local to the crack tip. It does appear though, if 
large scale yielding would take place, crack closure would have a 
lessening effect on the plate behavior. 

A further indication of greater constraint of yielding in the crack 
closure case is the relatively less crack tip blunting in the closure 
case than in the no closure case. The relative crack openings are shown 
in Figure (5- 21 )-. The craekb opening d,i s placements are; a. . - - 

normalized to the crack opening displacement at the plate center. It 
can be inferred from the data that for a given load, less plastic straining 
occurs in the closure case than in the no closure case. 
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FIGURE 5-20 

MOMENT REDISTRIBUTION FOR 
THE ELASTO- PLASTIC PROBLEMS 




v(r/a)/ v(l.O) AT 5 = 180 
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FIGURE 5-21 

CRACK TIP BLUNTING DUE 
TO EL ASTO- PLASTIC FLOW 
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In Figures (5-22), contour plots of w, — , and ~ are shown. 

ox oy 

Comparing these with the elastic results of Figures (5-2), (5-3), and 
(5-4) show that yielding has caused little change in behavior. This 
is due to the very local character of the crack tip yielding not affecting 
overall plate behavior. 
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! FIGURE 5-22 

CONTOUR PLOTS OF W, dw/dx a dw/dy 
AT THE MAXIMUM LOAD STEPS 
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CHAPTER VI 


FURTHER REMARKS 

A numerical solution to the problem of a centrally cracked plate 
subject to a circular bending field has been obtained. The distinguishing 
features of the solutions are that they allow isolation of the effects 
of crack closure and elasto-plastic material behavior. The numerical 
capability used to generate the solutions incorporates classic Kirchhoff 
plate theory assumptions with elasto-plastic work hardening material 
behavior. The results indicate that crack closure and elasto-plasticity 
are important considerations in the mechanics of through cracked plates 
subject to bending. 

6.1 Conclusions 

Specific conclusions with regard to this work are as follows: 

1. The effects of crack closure are significant in the mechanics 
of through crack bending and must be given consideration 

in the physical assessment of the problem. 

2. Neglect of the crack closure phenomena results in non- 
conservative calculations for crack opening displacements and 
bending stress intensities based on linear elastic fracture 
mechanics concepts. 

3. The most sensitive parameter to the crack closure phenomena is 
the transverse slope with respect to the direction normal to 
the crack face. (aw/ ay) 



92 


4. Neutral axis shifts, based on displacement measurements, 
clearly monitor the crack closure phenomena. As plastic 
yielding builds up, neutral axis shifts become less pronounced. 

5. Superimposing tension on exterior boundaries to preclude 
contact on the compressive^ surface due to bending does not 
model the crack closure phenomena in a physically realistic 
fashion. 

6. The Dugdale or strip model as applied to plate bending does 
not adequately model the elasto-plastic problem as it does 
not include the elastic core consideration. 

7. Plastic straining first occurs on the tensile surface of the 
plate and proceeds part-way through the thickness before 
yielding occurs on the compressive edge for the closure case. 

8. Inclusion of crack closure in the analysis inhibits plastic 
straining when compared to analysis neglecting closure, at 
least for the circular bending problem. 

6.2 Recommendations for Further Research 

The principal difficulty encountered in this work was the free 
surface approximation required by Kirchhoff plate theory. Knowles and. 
Wang [3] have shown that reformulation of the problem in terms of Reissner 
sixth order theory precludes this difficulty and results in solutions more 
compatible with the extensional case. Still, the crack closure phenomena 
has not been included and, based on these findings using Kirchhoff theory. 
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its significance needs to be evaluated using the higher order theory. 

A possible approach in this direction would be to employ Reissner plate 
theory in conjunction with the crack closure treatment used here. 

The inclusion of the geometric non-linearities with the associated 
coupling of in-plane and bending behavior needs to be evaluated especially 
with reference to the combined extension/bending problem. Formulation 
of this type of complex problem can be made more tractable through 
numerical determination of stationary points of functionals as used in 
this work. The conjugate gradient method may be formulated to find 
stationary values of non-linear problems as in [30] and [33]. 

Short of solving the complete three dimensional contact problem, 
experimental investigations are required to verify the assumptions made 
here regarding crack closure and to further evaluate the three dimensional 
effects of the crack closure phenomena. Particular interest lies in the 
evaluation of crack face warping, generation of through thickness stress 
variations, plate thickness to crack length ratio effects, and the 
influence of plate exterior boundaries on the stress and strain field 
data. 

Use of the solution capability to solve other bending problems of 
technical interest is also indicated. Elasto-plastic solutions to plates 
with patterned cutouts and under various boundary conditions are possible. 
Using the linear constraint equation concept, connection matrices between 
structures can be developed that model many different boundary flexibilities. 
This is a useful capability both in the practical engineering sense and as 
a research tool . 
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APPENDIX A 

CONSTITUTIVE EQUATION FOR PLANE STRESS 
For plane stress, the constitutive relation of Equation (3-9) is: 
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APPENDIX B 

BICUBIC HERMITE INTERPOLATION FORMULAS 

Displacement rate vector { u } and bicubic Hermite interpolation 
formulas used in Section IV. 



• 9 U 

u . . ■= — ; H . . = Hermite Interpolation Formula 

3X J at corner i,j 1J 

Displacement rates { v} and {w} are similar to {u}. 
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APPENDIX C 

ON THE DERIVATION OF FIELD EQUATIONS CHARACTERIZING 
THE ELASTO-PLASTIC BEHAVIOR OF VARIOUS SITUATIONS 

Using Theorem I developed in Chapter III, it is possible to 
derive field equations of el asto-pl astic flow for various physical 
situations. Using the strain-rate, displacement-rate relation 

‘ij ■ %( “ij + (c - 1) 

Equation (3-18), which is the condition for equilibrium, can be written as 

v 

J E ijkn“i ,j s "k,« dV 'I V“k ds “ 0 (c ' 2 > 

Integrating Equation (c-2) by parts produces surface and volume 
integrals plus constants of integration. Using the fundamental lemma of 
calculus of variations, this procedure will lead to field equations and 
complete boundary conditions. 

In the following sections, theories for plane stress/strain, St. 

Venant torsion, simple beam theory, and plate bending theory are presented 
to show both the adequacy of Equation (c-2) and to provide a convenient 
reference for the various elasto-plastic field theories. 
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C-1 Plane Stress/Strain 

Referring to Figure C-1, the usual plane stress assumptions are 


a - a = a =0 
xz yz ZZ 
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The constitutive equations are of the form 
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where [E] is given in detail for plane stress and plane strain in 
Appendix A. The volume integral can then be written per unit thickness as 
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FIGURE C-l 

PLANE STRESS GEOMETRY FOR 
THE STRETCHING THEORY 
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Integrating by parts and applying the fundamental lemma of 
variational calculus results in the field equations 


E ll°xx + 2E 13“xy + (E 33 + E 12>Vy + E 13 4 xx + E 23*yy 


+ (E 11 v + E 13 >“x + (E 12 + E 23 > v y 


+ (E ]3 + E33 )(u y + v x ) - 0 

.x ,y J 


(c-7) 1 


and 


E 22^yy + 2E 23^xy + E 33^xx + E 23**yy + ( E 33 + E 12^xy + E 13 u 


xy 


xy ^“xx 


* < E 22 „ + E 23 >*y + < E 12 „ + E 13 >x 


,X 


,X 


+ (E 33 + E 23 )(; x + V = 0 

>x ,y 


(c-8) 1 


with the complete boundary conditions at y = y Q or y - 0 
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* See J. L. Swedlow, "Character of the Equations of Elasto-Plastic Flow 
in Three Independent Variables," Int . J. Non-Linear Mechanics , 

Vol . 3, pp. 325-336, 1968, for verification of these equations. 
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02 St. Venant Torsion 

The usual St. Venant torsion assumptions using the notation of 
Figure 02 are 
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and kinematically 
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The constitutive relation is 
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TORSION THEORY 
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This can be inverted to be in the form 
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The strain-rate, displacement-rate equations are 




Equation (c-2) can now be written as 
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Integrating Equation (c-16) and (c-17) by parts and applying 
lemma results in the following theory for St. Venant torsion 
plastic rod; 

Differential equation 
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subject to the boundary conditions at z = z. 
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* See J. L. Swedlow, "Character of the Equations of Elasto-Plastic Flow 
in Three Independent Variables," Int . J_. Non-Linear Mechanics , 

Vol . 3, pp. 325-336, 1968, for verification of these equations. 
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C-3 Simple Beam Theory 

The usual assumptions accompanying simple beam theory consistent 
the notation of Figure C-3 are: 

1. Planes remain plane; i.e., the total strain distribution 
through thickness of the beam is linear. 

2. Effects of z-direction are neglected. 

3. Transverse deflection is a function of x only. . 

From assumption 1, the strain, displacement-rate equation can be 
written as 
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where u implies the extensional displacement associated with the mid 
plane of the beam. 

The constitutive relation of the plane stress problem reduces to 


The volume integral of Equation 
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GEOMETRY AND LOADING FOR BEAM THEORY 
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Defining 


y 0 

2t J E n dy 

-y„ 


y Q yo 

C 2 - 2t J E^ydy C 3 = 2t J Ell y 2 dy 


-y. 


-y. 
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2 2 
Note that for elasticity, E-|, = E/l-v a constant. Hence, C-| = 4Ety Q /l-v , 

C 2 = 0, and C 3 = EI/(l-v 2 ). Since C 2 = 0, the bending and stretching 
modes decouple in the elastic case. However, the coupling between in- 
plane and transverse displacements is present in elasto-plastic beam 
theory. 

The surface integral is written as 
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Integrating Equations (c-24) and (c-25) by parts twice and applying 
the fundamental lemma, the resulting differential equations theory for 
elasto-plastic beams is 
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subject to the boundary conditions at x 

= L 




r 

C ! a x- C 2*xx- F L ' 0 

or 

.S 

u 

Speci fied 

< 

(C 2°x>x - < C 3*xx>x + 1 L ■ 0 

or 


Specified 


, - (c 2 “x) + c 3\x ; \ - 0 

or 


Specified 


(c-27 ) 


at x = 0 




.S 

- C..U + C 0 v + F - 0 

lx 2 xx o 


.S 


or u S Specified 


- (C ? u ) + (C 0 v ) - Q„ = 0 or v Speicified 

^ X 


3 xx 


• S 


C 9 u - C.v + M =0 
^ x 3 xx o 


or v Specified 

A 


(c-28) 


For elastic problems as noted, = 0 and the above equations 
reduce to uncoupled elastic simple beam and uniaxial tension problem. 
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C-4 Plate Bending Theory 

Elasto-plastic plate bending theory may be derived by making 
kinematic assumptions similar to Kirchhoff bending theory. These 
assumptions include: 

1. Deflections are small compared to the thickness of the plate 
and the slope is small compared to one. 

2. State of plane stress exists in the plane of the plate, i.e.. 



3. Total strain distribution is linear through the plate thickness. 
For elasto-plastic plate theory, it is necessary to include the in- 
plane or membrane behavior as they do not decouple from the bending 
behavior as in elasticity. The plate and coordinate system are shown in 
Figure C-4. 

The strain displacement relations are 
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x x xx 
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y y yy 
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and the constitutive law is identical to the plane stress case given in 
Equation (c-4). 
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GEOMETRY AND LOADING FOR BENDING- 
STRETCHING THEORY 
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Using Equations (c-29) and (c-4), the volume integral of Equation 
(c-2) can be written as 
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+ 2E 23^y • ^^“y + 4 x ' 2Z \y ]2 + E 3 3 t V V 


(c-30) 


Since Ej. is the only quantity that is a function of z, the inte- 
gration can be reduced by defining 


E.. = T E. .dZ C.. * r ZE,.dZ D. . ° [ Z 2 E, .dZ 

ij J ij ij J ij ij J ij 


(c-31) 


-z 


-z 


With E.j, C^. , and D_ introduced into Equation (c-30) and expanding 
terms, the volume integral may be written as 

y x 

o o 

l E ijk^i.o 4C, k,x dv = I J ( cr i l °X + r i2*y + F 13 { “y + V ]sCl x 


o o 


+ ^Vx + E 22V + E 23<“y + \ )W y + [E l 3 “x + E 23*y + E 33 ( V "x^'V V 


[C n w + C 10 w + 2C,,w ]fiu v - [C,,w + C ft + 2C 00 ft ]6v 

11 xx 12 yy 13 xy x 12 xx 22 yy 23 xy y 


[ C 13"xx + C 23*yy + 2C 33*xy ]S(0 y + *x> V C 12V C 13 ( V V ]5 “ X x 


(This equation is continued on the next page.) 
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[C 12 U X + C 22 V y + C 23 (u y + v x )]i "yy ' 2[C 13 U x + C 23% + C 33 ( “y + 4 x^ Sw «y 


+ [D ll\x + D 1 2^yy + 4D 13V s *xx + tD 12«xx + D 22"yy + 4D 23 ,:, xy ;|6 "yy 


+ 4 t D l 3 *xx + D 23*yy + ^xy 34 V dXdy 


(c-32) 


The surface integral of Equation (c-2 ) is 

v y ° 

f T.du.dS = + f [N 6u + N dv + M 6w + M dw - Q dw]dy 

•lii J x xy xx xy y x 


f [N dV + N du + M dw + M dw - Q 6w]dx 
J y xy 


y y yx x y 


x o ^o 


j* J pdwdxdy 


o o 


where 


(c-33) 


| N y > 
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r 'a 
a 

x 


xy 


* I 


M 


-z 


J 


)dz, / 




V 


O 

J 2 


-z 




Integrating by parts and applying the fundamental lemma of varia- 
tional calculus results in the field equations which couple the in- 
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P aCeTO " ts the c . s r „ 

f J r,,ese nations 

HVx * 2f0, 2 * 2t J; 

33 x*yy Vyyyy + Sfn i + . 

23 Wy + 0 I3 V y J 

* T&) 1 + n 

I I U l 9 + 40 T • 

.xx 12 vi/ 13 Jw + rn 

•w ,3 ,xy *x 10] 2 + 0 +4D 

.XX 22 V1/ ^ u 2 3 Jw 

,yy ,xy yy 


+ 4 C°„ *ft 


13 + ^?o +0 v 

.yy 33 , xy J "xy ' c n%< x - 3 C ) 3 .u, 

“'"‘■•Vvi-.. 


xyy yxx 3 " ^23^v, - C v 

yyy 3 v 


- 3 C 9o v 


XXX 


fc-35 


'23*xyy • C 2oV wl . r c 

3 2 yyy K„^ tc 3 . 

,Xy < 


~ +JC . 

’X* ?3 ,yy X, xy 1[u y +\] 

£ n\ x + an ij + f • 

^ £ 33 “yy M £ 33 + F , 2 )v + f . 

3 xy E,3V XX ♦ F 23 v yy 


+ f£ n +F r 

13 A * [E 


Ju + rr _ 

,y X t £ , 2 +£ JV + TF - 

»x v y ‘•^70 +£ t / • 

,y 3 ,x 33 3 ‘V v ) 

.a ,y y X 7 


c "%x - 3 C, 3 w . . 

3 Xxy (C,J +2C ^ 

3 xyy % Vy 


(c-36J 



C-18 


E 22*yy + 2E 23*xy + E 33°xx + E 23“yy + [E 33 + E 12 ]Ct xy + E 13 6 xx 


+ CE 22 v + E 23 A + [E 12 +E 13 A + [E 33 X +E 23 A A 3 

9J 9J 9* 9J 


' C 13 w xxx ‘ ^ C 12 + 2C 33^ w xxy " 3C 23 w xyy + C 22*yyy 


[C-j 2 + C-J3 ]& xx - + 3^ 


,x 


23 y 33 ,x ** 


[C 22 + C 23 ] \y ' 0 

9J 9* 


(c-37 ) 


The complete boundary conditions are: 

At the corner point x = x Q and y * y Q 

M xy + M yx + 4 ^1 3^xx + °23^yy + °33^xy^ ' 2 C C 13 fl x + C 23V C 33A + °x )] = 0 


or w Specified 


(c-38) 


Along the line x = x q for all y, there are four conditions as follows: 
*11 °x + E 12*y + E 13 (u y + V " C ll W xx “ C 12 w yy " 2C 13 W xy ' 2Z^ = 0 


or u Specified 


(c-39) 
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_ , . N 

r i3“x + f 23*y + E 33<“y +ii X> ‘ C 13*xx ' C 23"yy ' 2C 33*xy ' # ' 



2 * C 33 V x^ 


or w Specified 


D ll*xx + D 12*yy + 4D 13*xy " [C llV C 12 v y + C 13^y + v x^ = M x 

or w Specified 
X 


Along the line y s y Q for all x, 

f 13“ x + f 23 ; y + r 33 (0 y + V - [C 13“xx tC 23*yy + 2C 33V ' 2^ = 

or 0 Specified 


r !2°x + Vy + £ 23<V } x> - [c l2*xx + c 22*yy + 2C 23V " 2% 


0 

(c-40) 


(c-41 ) 


(c-42) 


0 

(c-43) 

0 


or v Specified 


(c-44) 



a* [4(D 13 w j + . 

xx * h ' u 23 \ v ) ■ + 


4 °23%y * [(C 12 UJ,( t 
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y '-32V J 


y * + (c & V„ + «a* x j 


+ 2(C,,i 


13Ux) x + + 2fC„(u + v 


33 { V vJJ 


j' V J X + * 0 


or 


w Specified 


" xx ♦ V yy * 


[C '^ c 22V^iu y+i)] 
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(c-45) 


or 


w y Specifi ed 


(c-46) 
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APPENDIX D 

BOUNDARY CONDITIONS AT CORNER POINTS IN KIRCHHOFF PLATE THEORY 


As noted in Equation (c-38), there are reactive forces present at 
corner points in fourth order plate theory. In this appendix, the 
boundary conditions at points where boundaries are discontinuous are 
derived via the variational principal stated in Equation (c-2). The 
plate under examination (Figure D-l) has no applied in-plane loadings and 
is assumed to be elastic. 

The volume integral of Equation (c-32) is reduced to include only the 

transverse deflection noting that for elastic behavior* = 0. However, 

r v 

the surface integral, T.u.dS, which includes the bending moments M 

l l n 

and M nt , and the transverse shear force Q n along with the transverse 
pressure p, becomes 

J Tu.dS = - j 

S a 

By setting the first variation of the volume plus surface integrals 

to zero, the boundary conditions for the case where surface tractions 

are applied over the surface Sa and Sa in Figure D-l can be found. At 

1 2 

the intersection of Sa and Sa , either the transverse deflection w is 

1 2 

zero, or 

[(w xx - w yy )(sin 2 m-, -sin 2« 2 ) 

+ M ntJ ' 

on S-| on S 2 


2w (cos 
xy 


2a j -COS 




M nt 


M 

n 3n 


M Di- 
nt at 


Q n w 


dS 


(d-l) 


0 (d-2) 



FIGURE D-l 

PLATE GEOMETRY FOR CORNER EXAMINATION 



NOTES : 

1. APPLIED TRACTIONS ON S^, a S a2 

2. APPLIED DISPLACEMENTS ON S u 

3. n a t REFER TO NORMAL AND TANGENTIAL 
DIRECTIONS AT A POINT 

4. X 8 Y ARE CARTESIAN COORDINATES 
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The various corner configurations of Figure D-2 are examined in 
light of Equation c-47 and are listed in Table D-l . It is noted that 
the reactive conditions at the intersection vary from zero for a smooth 
intersection to a maximum for the exterior corner, back to zero for the 
crack tip, and to a minimum for the interior corner. 
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TABLE D-l 

Corner Reactions in Elastic Plate Bending Theory 


FIGURE 

“1 

°2 

M . 
nt 

on 

s °i 

M nt 
on 
So 2 

CORNER REACTION 
EQUATION (d-2) 

Smooth Intersection 
Figure A-6.a 
D-2 • 

a 

a 

M nt 

M nt 

' 0 

Exterior Corner 
Figure A-6.b 
D-2 

0 

90° 

M xy 

M yx 

M - M = - 2D ( 1 -v)w 

yx xy v J xy 

Crack Tip 
Figure A-6.c 
D-2 

90° 

-90° 

• 

M 

yx 

M yx 

0 

Interior Corner 
Figure A-6.d 
D-2 

90° 

0° 

M 

yx 


M - M = - 2D(1 -v)w 

xy yx xy 
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